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Abstract 

The influence of turbulent effects on a fluid flow through a (pseudo) 
porous media is studied by numerically solving the set of Reynolds-averaged 
Navier-Stokes equations with the K-e model for turbulence. The spatial 
domains are two-dimensional rectangular grids with different porosities 
obtained by the random placing of rigid obstacles. The objective of the 
simulations is to access the behavior of the generalized friction factor 
with varying Reynolds number. A good agreement with the Forchheimer's 
equation is observed. The flow distribution at both low and high Reynolds 
conditions is also analyzed. 

1 Introduction 

The flux of fluid through microscopically disordered and macroscopically ho- 
mogeneous porous media is a function of the pressure gradient, which for low 
Reynolds conditions is well approximated by a linear relation known as Darcy's 
law §, 

VP = -f, (1) 

where P and u are the averages of pressure and velocity, respectively, over cross 
sections perpendicular to the fluid flow, /i is the fluid viscosity and k is the the 
global permeability. This was formulated by the French engineer Henri Darcy 
based on experimental observations of the pressure driven flux of water through 
sand filters. However, experimental observations and numerical simulations have 
shown that the relation between pressure gradient and mean velocity is indeed 
nonlinear, and can be well fitted by the so-called Forchheimer's equation J|, |j. f| 

- VP = afiu + f3pu 2 , (2) 

where the coefficient a corresponds to the inverse of the permeability and (3 is 
usually called the inertial parameter. The transition from the linear (Darcy's 
law) to the nonlinear regime occurs gradually as the Reynolds number in- 
creases. In order to understand the connections between the pore morphology 
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and the macroscopic behavior of the system it is necessary to analyse the spatial 
pore distribution and relate it to the mechanisms of momentum transfer (iner- 
tial, viscous and turbulent). Previous studies have been conducted in this direc- 
tion leading to successful prediction of permeability coefficients of real porous 
materials ||.The role of inertial effects over such a transition has been analyzed 
through computational simulations in |^|, |(| within laminar flow regime through 
the pore space. 

The random aspect of the pore distribution induces a highly heterogeneous 
local flow which becomes turbulent at high Reynolds regimes. Following the 
settings of previous numerical works J5], |(| we are going to analyse the combined 
effects of inertia and turbulence in the behavior of flux of fluid through two- 
dimensional pseudo-porous systems, whose pore connectivity is based on the 
ideas of site percolation disorder. To set the geometry for the simulations we 
start with a square lattice of 64 x 64 obstacles from which cells are removed at 
random until we get a pseudo-porous medium with a prescribed void fraction. 
To minimize the end effects over the flow field in the pore region we attach a 
header region in the inlet and a recovery region in the outlet, both of them free 
of obstacles. The fluid is assumed to be incompressible and Newtonian and we 
adopt the n-e model of turbulence. The velocity field at the inlet is assigned 
a constant value normal to the boundary. At the outlet we fix a null gradient 
velocity boundary condition. 

One goal of our simulations is to analyse the exactness of Forchheimer's 
equation eq.(^|) in representing the behavior of turbulent fluid flow in porous 
media. We rewrite eq. (||) in a form that has been successfully used to correlated 
experimental data from a large variety of porous materials and flow conditions, 
viz., 

/-I-5?. (3) 

where / = —VP/((3pu 2 ) and Re' = (3pu/(afi). The flow localization over the 
pore space is quantified by the partition function [[| || , 

(4) 

where n is the total number of fluid cells and qi — Ei/ Y^j=i Ej, with Ei being 
the kinetic energy associated with the fluid cell i. This is a measure of the 
uniformity of the flow through the available pore space. If the flow is evenly 
distributed over the entire number of cells, tt attains the value 1. On the other 
hand, if the kinetic energy is concentrated over few cells the function n would 
be close to the value l/n. 

2 Turbulent Flow and the tz-e model 

To account for the turbulence we use the Reynolds averaging approach of de- 
composing the variables into a mean (ensemble-averaged or time-averaged) and 
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a fluctuating component 

<f>i = <j>i + 4>[ (5) 
where (j) represents a generic variable. For the velocity components we have 

Ui = Ui + u-. (6) 

Substituting the decomposed expressions for the velocity and the pressure into 
the Navier-Stokes equations for incompressible fluid and dropping the overbar 
on the mean velocity and mean pressure we get the Reynolds-averaged Navier- 
Stokes equations 
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In order to close these equations, the Reynolds stress term —pu^u'j must be 
modeled, and a common method employs the Boussinesq hypothesis to relate it 
to the mean velocity gradients: 



- PUiU'j = IM 
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(9) 



In the case of the K-e models, two additional transport equations (for the 
turbulent kinetic energy, k, and the turbulent dissipation rate, e) are solved, 
and fi t is computed as a function of k and e, 



(10) 



with C M = 0.09 a constant. The equations for k and e are 
Dk 
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represents the generation of turbulent kinetic energy by the mean velocity gra- 
dients. To evaluate Gk in a manner consistent with the Boussinesq hypothesis, 
the Fluent code is implemented with 



G K = fi t S 2 
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where S is the modulus of the mean rate-of-strain tensor, defined as 



S — y 2Sij S, 



with the mean strain rate Sij given by 



1 / diii duj 

2 \dxj dxi 



) 



We use the Fluent default values for the constants appearing in the above 
equations, i.e., C% e = 1.44, Ci^ = 1.92, a K = 1.0 and cr £ = 1.3. 

We use the software Fluent to solve these equations in our simulations, which 
has implemented the n-e model of turbulence. 

The criteria for convergence we used in the simulations is defined in terms 
of the residuals which provide a measure of the degree to which each of the 
conservation equations are satisfied throughout the flow field. Residuals are 
computed by summing the imbalance in each equation for all cells in the domain. 
The residuals for each flow variable (e.g. velocity, pressure, etc.) give a measure 
of the error magnitude in the solution at each iteration. In general, a solution 
can be considered well converged if the normalized residuals are on the order of 
10~ 3 . In all of our simulations, convergence is considered to be achieved only 
when each of the normalized residuals fall below 10 -5 . 

3 Experiments 

We performed experiments with three values of porosity (e = 0.7, 0.8, 0.9) which 
are greater than the critical percolation porosity. For each value of e, we gener- 
ated ten different grids as described above, and with each one we simulated the 
flows using 17 values of viscosity 

The velocity and pressure fields in the pore, header and recovery regions 
were numerically obtained through discretization by means of the volume finite- 
difference technique^. In each simulation, after we reached a converged flow, 
we did a post processing to compute the difference between the averaged values 
of the pressure on the cross sections at the beginning and at the end of the pore 
region. We computed also the participation function n. Then we computed the 
average over the ten realizations and construct the curves of fix AP and fiXiz. 
Using a linear regression on the graph of fix AP we determined the constants 
a and /3 and then computed the values of Re' and /. The results are presented 
through the graphs of Re' x (/ — 1) and Re' x it. 

4 Results and conclusion 

Our simulations may be summarized at two graphs, one for Re' x (/ — 1) and 
the other for Re' x tt. At both graphs we can see an agreement with previous 
simulations |5|, |6). In the graph of Re' x (/ — 1) we can see that the Forchheimer's 
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equation (dashed straight line) is indeed a good fit for the numerical results up 
to Re' w l,for e = 0.7 and e = 0.8, and up to Re' w 10 1 for e = 0.9. 

The inclusion of turbulence makes it possible to extend the simulations to a 
wider range of Reynolds number, including values that would possible present 
convergence problem if considering the laminar Navier-Stokes equations. This 
is due to the addition of the turbulent viscosity which increases the energy 
dissipation. 

The higher the porosity the bigger is the participation function, which ex- 
presses the fact that the flow gets more uniform as the porosity increases. Also, 
we can see that the variation of it increases as the porosity gets higher. We 
can also identify two plateaus: one for low and the other for high values of Re'. 
This is an indication that for low velocities (low Re') conditions, the flux is 
mainly through a few preferential channels, the number of which increases with 
the porosity. As the velocity (Re') increases more channels get to be used in- 
creasing the participation function n, which however, has an upper limit smaller 
than unity. 
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Figure 1: Plots in log-log scale showing the relation between the generalized 
friction factor / — 1 and the modified Reynolds number Re' for the experi- 
ments with porosities e = 0.7, 0.8 and 0.9. The dashed straight line represents 
Forchheimer's equation (eq.(Q)). 




Figure 2: Log- linear plots of the modified Reynolds number versus the partic- 
ipation function n for the porosity values e = 0.7, 0.8. 
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